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Abstract 

In this paper, we consider the numerical solution of some nonlinear poroelasticity problems that are 
of Biot type and develop a general algorithm for solving nonlinear coupled systems. We discuss the dif¬ 
ficulties associated with flow and mechanics in heterogenous media with nonlinear coupling. The central 
issue being how to handle the nonlinearities and the multiscale scale nature of the media. To com¬ 
pute an efficient numerical solution we develop and implement a Generalized Multiscale Finite Element 
Method (GMsFEM) that solves nonlinear problems on a coarse grid by constructing local multiscale 
basis functions and treating part of the nonlinearity locally as a parametric value. After linearization 
with a Picard Iteration, the procedure begins with construction of multiscale bases for both displacement 
and pressure in each coarse block by treating the staggered nonlinearity as a parametric value. Using 
a snapshot space and local spectral problems, we construct an offline basis of reduced dimension. From 
here an online, parametric dependent, space is constructed. Finally, after multiplying by a multiscale 
partitions of unity, the multiscale basis is constructed and the coarse grid problem then can be solved for 
arbitrary forcing and boundary conditions. We implement this algorithm on a geometry with a linear 
and nonlinear pressure dependent permeability field and compute error between the multiscale solution 
with the fine-scale solutions. 


1 Introduction 

The applications of mechanics and flow in porous media are wide ranging, as are the challenges involved in 
simulating some of these problems in nonlinear multiscale contexts. This is particularly true in geomechan¬ 
ical modeling where relevant phenomena may be highly nonlinear, for example in the setting of enhanced 
production and environmental safety concerns due to overburden subsidence and compaction [T^ |2Q] . An¬ 
other of the central challenges is the multiscale nature of the media considered in geomechanics problems. 
Heterogeneity of rock properties should be accurately accounted in the geomechanical model, and this re¬ 
quires a computationally costly a high resolution solve. Moreover, due to the multi-physics nature of the 
problems, they may involve highly nonlinear relations. This then makes the further requirement of many 
iterations in a Newton or Picard linearization. Thus, we propose a multiscale method to attempt overcome 
some of these challenges. The central idea is to linearize in a Picard iteration, and treat the nonlinearities 
as a parametric value as utilized in m and references therein. 

As noted in [S], the basic mathematical structure of the poroelasticity models are usually coupled equa¬ 
tions for pressure and displacements known as Biot type models |24j . The pressure equations are a parabolic 
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equation coupled to a time derivative of volumetric strain. While the mechanics equations are are given by 
quasi-static elasticity equations and is coupled by gradients of pressure. In this work however, we focus on 
the possible nonlinear couplings of the Biot model. There are a myriad of physical and modeling reasons 
to add nonlinearity to the Biot equations, however, we will primarily focus when the permeability field and 
elasticity tensors depend nonlinearly on pressure and displacements and their gradients. This is due pri¬ 
marily to us wanting to focus on the nonlinearities effects on our GMsFEM, as nonlinearities in lower order 
derivative will not interfere with the construction of the local multiscale basis functions. 

Nonlinear Poroelastic models of this type have been explored in the literature to incorporate higher order 
physics considerations. For example, when the viscosity of the fluid heavily depends on the fluid pressure 
we may obtain relations of permeability of the form 


K{x;p) 


k{x) 

fj.{x;p)' 


Here k is the absolute permeability, and p{x;p) is the pressure and spatially dependent viscosity. This can 
occur when their are very high pressure gradients |29) . In the setting of complex geomechanical interactions 
used a relationship between permeability and volumetric strain of the form 


iF(a;; V ■ u) = A{x) exp(H(a;)V • u), 


where A,B are determined constants and V • u is the volumetric strain. Further in |15j . the porosity 
(j) also depends linearly on (Vp, V • u), however, this is multiplied throughout generating a nonlinearity. 
In the context of fractured reservoirs, permeability is often computed via the so called ’’cubic-law” through 
channels and this may be coupled in orientation and magnitude via the displacements in a nonlinear way [30] . 
With this GMsFEM, we propose a method to efficiently compute solutions to these nonlinear poroelasticity 
problems with the heterogeneous multiscale properties. 

As noted in the prequel [9] , there are many very effective multiscale frameworks that have been developed 
in recent years. There are rigorous approaches based on homogenization of partial differential equations mM- 
However, these approaches may have limited computational use. Examples of computational approaches in¬ 
clude the Heterogeneous Multiscale Method (HMM) [5^ , approaches based on the Variational Multiscale 

Method (see m), where coarse-grid quasi-interpolation operators are used to build an orthogonal splitting 
into a multiscale space and a fine-scale space [28]. In this paper, we will use the Generalized Multiscale Finite 
Element Method and its extension to nonlinear poroelasticity problems in the framework of |10j . Specifically, 
to handle multiscale nonlinear problems, we combine ideas of model reduction, whereby the nonlinearity is 
replaced locally by a parameter space and offline and online spaces are generated. For a broad presentation 
of these methods we refer the reader to [T]. 

The paper is organized similarly to [9|, as follows. In Section 2 we provide the mathematical background 
of the nonlinear poroelasticity problem. We will introduce the Biot type model and highlight where the 
heterogeneities primarily occur. We again note that the nonlinearities in our model are in the permeability 
and elasticity tensor as these are second order derivative terms. In Section 3, to outline the difficulties 
in full direct numerical simulation we introduce the fine-scale discretizations using coupled time-stepping 
schemes and a Picard iteration technique for linearization. In Section 4, we present our nonlinear GMsFEM 
algorithm and outline its construction procedure. Finally, numerical implementations are presented in Section 
5. Using the GMsFEM, we compare the multiscale solution to fine-scale solutions and give error estimates. 
We will present two different examples with permeability being linear and nonlinear with respect to pressure. 
Additionally, we will implement and discuss different snapshot spaces and coarse-grids choices, and its relation 
to enrichment and the error. 
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2 Problem formulation 

We denote our computational domain n C to be a bounded Lipschitz region. We consider a general 
nonlinear poroelastic system where we wish to find a pressure p and displacements u satisfying 

— div {C{x;u,\/u,p,\/p)e{u)) — a\/p = 0 in fl, (la) 

-div{K{x;u,\7u,p,\7p)\7p)+a^^^^ + = / in 11, (lb) 

with initial condition for pressure p{x,0) = pq. We write the boundary of the domain into four sections 
i9n = riUr 2 = r 3 Ur 4 . We suppose the following boundary conditions on each portion 

C(x\u,Vu,p,'S/p)e{u) ■ n = a; € Ti, u = ui, a; € r2, 


and 


dp 

-K{x;u,\/u,p,\/p)—=0, a: e Ta, p = pi, a; S r 4 . 


Here the symmetric strain is written as e{u) = | (Vu + Vu^) and we write 

C{x] u, 'S/u,p, 'S/p)£{u) := C(x; it, Vu,p, 'S/p) : £{u) 


to mean the double contraction of a 4-tensor with a 2-tensor. 

As in the linear case, the primary sources of the heterogeneities in the physical properties arise from 
C{x]u,S/u,p,S/p), the elastic tensor, and K{x;u,S/u,p,S/p), the permeability. In this setting, we suppose 
these heterogenous parameters can depend in p and u and their gradients in complicated nonlinear ways. 
Further, we will denote M to be the Biot modulus, v is the fluid viscosity, and a is the Biot-Willis fluid-solid 
coupling coefficient. Here, / is a source term representing injection or production processes and n is the unit 
normal to the boundary. Body forces, such as gravity, are neglected. 

Remark; Note that one could also add nonlinearities in the coefficients a and M, however, these 
correspond to lower order terms with respect to derivatives. Therefore, these will not contribute to the local 
problems in the GMsFEM. Hence, we will consider them to be constant throughout. 

We recall the setting when these relations become linear. In the case of a linear elastic stress-strain 
constitutive relation we have that the stress tensor and symmetric strain gradient may be expressed as 

Ce{u) = 2fi£{u) + A div(u) I, 


where p, A are Lame coefficients, I is the identity tensor. Note here this p is not to be confused with what 
is often used as a parameter. Above we have absorbed into the nonlinear permeability coefficient the fluid 
viscosity and in the case of linear permeability, we have 



k being absolute permeability. 

The nonlinear poroelasticity problem Q, can be written in operator matrix form: 

A{u,p) + aGp = 0, (2) 

^{Sp + aDu) + B{u,p) = f, (3) 

where 

A{u,p) = — div {C(x\ u,S/u,p, Vp)£(u )), B{u,p) = — div {K{x; u, S/u,p, Vp)Vp), 
and G and D are gradient and divergence operators and S = 
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3 Fine-Scale Discretization 


We will now present fine-scale approximation and nonlinear solution methods for the above system. We will 
motivate the need for a multiscale method due to the nonlinearity and the heterogeneity of the poroelasticty 
problem. To approximate the solution to Q on fine-scale grid we will utilize a standard finite element 
method. The corresponding nonlinear variational form of the continuous problem written as 

a{u, p, v) + g{p, v) = 0, for all v G V, (4) 

^ =(/>'?)> forallgSQ. (5) 

for u G V, p G Q where 

V = {v G : v{x) = ui,x G r 2 }, Q = {q G : q{x) = pi,x G r 4 }, 

and the test spaces with homogeneous boundary conditions are given by 

V = {v G : v{x) =0,x G r 2 }, Q = {q G : q{x) =0,x G r 4 }. 

We define the following nonlinear forms 

a{u,p,v) = / {C{x;u,\7u,p,\7p)e{u),e{v))dx, (6a) 

Jn 

b{u,P,q)= [ {K{x;u,Vu,p,Vp)Vp,\7q)dx, (6b) 

Jq 


and bilinear and linear forms 


c{p, q) 


1 

M 


pqdx, 


9{p,v) 


/ a{\/p,v)dx, 

Jn 


and 

d{u,q) = / adivuqdx, {f,q)= / fqdx. 

Jn Jn 

Here (•, •) under the integrand denotes the standard inner product. In Section]^ we will discretize the spaces 
using a fine-scale standard FEM and denote them 14, Qh and Vh,Qh, h being the fine-grid size. The FEM 
using these spaces will serve as a reference solution for our GMsFEM outlined in Section 

Nonlinear Solve: We will first consider the time discretizations of the above system, then will discuss 
resolving the nonlinearity. This discretization leads to several possible couplings between time-steps and the 
two equations of linear poroelasticity nmH]. However in the nonlinear case we will only consider the fully 
coupled scheme. We proceed by introducing for the nonlinear fully coupled time derivative operators and 
then the Picard iteration for the linearization of the nonlinear operators. 

The standard fully implicit finite-difference scheme, or coupled scheme, can be used for the time- 
discretization and is given by 


a{u^+\p^+\v)+g{p^+\v)=Q, (7a) 

+b{u-+\p'^+\q) = {f,q), (7b) 

with u" = u{x, tn), p" = p{x, tn), where = nr, n = 0,1,..., Mt, Mtt = T and r > 0. 

We will now consider nonlinear solve in space after time discretization by the fully coupled scheme Q. 
One could rewrite Q as a nonlinear system each time step and use a Newton solver, however, for our 
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GMsFEM we prefer to use a linearization based on Picard iteration. Indeed, we may linearize Q given 
{pj,Uj) from a previous iteration step we write 

a{uj+i,pj+i,v) K. a^{pj\Uj+i,v) := / {C{x]Uj,Vuj,pj,Vpj)e{uj+i),v)dx, 

JQ. 

h{uj+i,pj+i,q) Kih^{pf,Pj+i,v) := / {K{x]Uj,Vuj,Pj,Vpj)Vpj+i,Vq)dx, 

Jn 

where pj = {uj,Vuj,pj,Vpj). We choose this notation in part to emphasize this may be viewed as a 
parameter in offline phase of the GMsFEM. 

Fixing the time-step at (n -|- 1), and taking pj = as data from the previous 

iteration. For j = 0,1, 2,..., we wish to find such that 

{Po ; u'jtl . + aip^jXi > '^) = 0> (8a) 

+b^{ppp]tlq) = iLd), (8b) 

Once the desired convergence criteria is reached, we can set the terminal ,p^~^^) as previous time data. 
We then return to the algorithm time-stepping and continue the iterative linearization until the terminal 
time. Note that this process can also be used in an appropriate nonlinear generalization to a fixed stress 
splitting [T^fTS] . 


4 GMsFEM for nonlinear poroelasticity problem 

We will present the offline and online multiscale basis construction in the fluid or pressure solve then its 
construction in the mechanics or displacement calculation step in this nonlinearly coupled formulation. 
Similar to the presentation outlined in [9], however, we will focus on the effects of the nonlinearities on 
the method. Observing the linearized formulation (|^, we see that we may consider the nonlinearity as 
parametric values we are able to successful design a GMsFEM for this nonlinear problem. In this way, we 
are able to construct an online-offline multiscale basis with respect to this nonlinearity. We now outline the 
general procedure of the GMsFEM algorithm. 

We begin briefly with some standard notation. The overall fine-scale model equations will be solved on 
a fine-grid using spaces Vh, Qt and 14, Qh, and will be used for our reference solutions. We now introduce 
the coarse grid. Let be a standard conforming partition of the computational domain 17 into finite 
elements. The fine-grid, 7h can be taken as a refinement of the coarse-grid. We refer to this partition as the 
coarse-grid and assume that each coarse element is partitioned into a connected union of fine grid blocks. 
We use {xi]f^i, where N is the number of coarse nodes, to denote the vertices of the coarse mesh T^, and 
define the neighborhood of the node Xi by 

j 

See Figure[l]for an illustration of neighborhoods and elements subordinated to the coarse discretization. We 
emphasize that the use of uJi is to denote a coarse neighborhood, and we use K to denote a coarse element 
throughout the paper. 

For global coupling we use the linearized continuous Galerkin (GG) formulation to find {u^XXp]+i ) ^ 
( 1 ^ 0115 14n) such that 

(Mj ; UjXi . v) + g(PjXi > '^) = 0> (9a) 

+b^(pp,p;;XiXd) = (f,d), (9b) 
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Figure 1: Illustration of a coarse neighborhood and coarse element 

where Qon and V^n denote the online spaces. The online spaces are spanned by multiscale basis functions 
V'm m) and for n + 1 time step and j-th iteration, each of which is supported in Wi 

i,m i,k 

The indexes m, k represent the numbering of these multiscale basis functions for pressure and displacements, 
respectively. Here the parameter fi represents the nonlinear dependence as in ([^. Recall that we may take 
Mi = from the previous time-step and will treat these variables as parametric values on 

each coarse patch. However, for simplicity we will suppose that the dependence is only on (u,p) of this 
nonlinearity. 

Remark: Note, the derivative dependent problems, with nonlinear couplings of (Vu, Vp), may be han¬ 
dled. However, due to the oscillation in these quantities, these terms may not be well approximated by 
constants on the coarse-grid level. Thus, we would need to have a more enriched parameter space than is 
utilized here. 

We now discuss further how we handle the parametrized nonlinearities. We assume that u and p are 
bounded above and below, i.e. U G [Umin,Umax] and p G [pmin,Pmax], where {Umin,Umax) and {pminjPmax) 
are pre-defined constants. These may be guessed initially based on initial data or a-priori estimates. The 
intervals [urmn,Umax] and [pmm^Pmax] are divided into N equal regions: 

'^min ~ ^0 ^ ^ ■■■ ^ 1 ~ '^max') 

and 

Pmin — PO Pi ^ PN—1 ^ PN — Pmax- 

Clearly, if necessary these domains can be partitioned in different number of regions, but for simplicity we 
suppose they are equal in number. For the parameter pj we take average values of and p’J^^ in each 
coarse region uJi. For average of a function we will use the notation 

f = T^f fdx- 

|Wi| 
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More specifically, we use fxj to represent the dependence of the solution on . The multiscale 

basis functions will be computed for a selected number of the parameter values fXj, j = 0 ,N at the offline 
stage and we will compute multiscale basis functions for each new value of for each oji at the 

online stage. 

Boadly speaking, the GMsFEM algorithm consist of several steps: 

• Offline computations: 

1. Generate the coarse-grid, ■ 

2. Gonstruct the snapshot space, used to compute an offline space, by solving many local problems 
on the fine-grid. 

3. Gonstruct a small dimensional offline space by performing dimension reduction in the space of 
local snapshots. 

• Online computations: 

1. In each time step and nonlinear iteration for current value of fXj in each Wj, we compute multiscale 
basis functions and construct online space by performing dimension reduction in the offline space. 

2. Use small dimensional online space to find the solution of a coarse-grid problem for any force term 
and/or boundary condition. 

We construct multisclate basis functions for pressure and displacements separately. We begin by consid¬ 
ering the pressure solve, then, the displacement solve. 


4.1 Multiscale basis functions for pressure 

In the offline computation, we first construct a snapshot space (5“„ap- Construction of the snapshot space 
involves solving the local problem for various choices of input parameters and various boundary conditions. 
These local spatial fields are used then used construct the offline space and the space consists of fields defined 
on a fine grid. There are a few options available when constructing the snapshot space and we will proceed 
with the two most natural ways. 

Snapshot Space 1: First, we propose a snapshot space generated by harmonic extensions of b^. For 
simplicity, we will omit the index i when there is no ambiguity. We thus define such that 

= 0 in w, 

( 10 ) 

Here (x) are defined by (x) = Si^k, VI € where J/i(a;) denotes the fine-grid boundary node on duj. 

This is done for each fixed parameter Hj, j = 0,...,N. 

Snapshot Space 2: Alternatively, we may use local fine-scale space basis functions within a coarse 
region and construct local snapshots by solving the following eigenvalue problem with natural boundary 
conditions 




Where 


= / (.K{x,^ij)V4>i,V(j)j) dx, / K{x,Hj)4>i<t>jdx, 

Jn Ja 


(j>i are the standard fine-scale basis functions, and for each fixed parameter values /ij, j = 0,...,N. 
Let li be the number of functions in the snapshot space in the region w, and define 

Q“nap = span{V’/7'' : 0<j<N}, 
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for each coarse subdomain lo. We reorder the snapshot functions using a single index to create the matrix 


RP = 

^snap 






where Mgnap denotes the total number of functions to keep in the snapshot construction. 

To construct the offline space Qoff, we perform a dimension reduction of the space of snapshots by using an 
auxiliary spectral decomposition. More precisely, we solve the eigenvalue problem in the space of snapshots: 

where 


dx = (i?Lap)^S^Lap, 

J UJ 


Here 


N 


K{x) ='^tjK{x,^j,j), 


i=i 


is independent of /ij and tj are prescribed non-negative weights. The main objective is to use the offline 
space to accurately construct a set of multiscale basis functions for each fxj in the online stage. At the offline 
stage the bilinear forms are chosen to be parameter-independent, such that there is no need to reconstruct 
the offline space for each 


We then choose the smallest eigenvalues from Eq. (12) and form the corresponding eigenvectors in 
the space of snapshots by setting for fc = 1,..., N^g^, where are the coordinates 

of the vector We denote the span of this reduced space as Q'^g. 

At the online stage, for a given parameter value /i, multiscale basis functions are computed based on each 
local coarse region The associated online space is the small dimensional subspace of the offline 

space. In particular, we seek a subspace of the offline space that can approximate any element of the offline 
space in an appropriate sense. In the the online stage the bilinear forms are chosen to be parameter-dependent 
and we use following eigenvalue problem 


^on^on ^ 


(13) 


where 

[ {K{x,y:)^<j^f,^<j>f)dx=[RlgfBRlg, 

J UJ 

M™= [ K{x,fi)cl,fr/dx = {RPgfMR^^g. 

J (jj 

Here B and M are the fine scale matrices corresponding to the stiffness and mass matrices for given p, and 


7?P — 


/ off / off 

Vl ; • ■ ■ ; 

off 


Finally, we multiply the partition of unity functions Xi by the eigenfunctions in the online space to 
construct the resulting basis functions 


fo'' 1 < * < and 1 < fc < , (14) 

where i/'™ = Y^j=i Xi is the standard linear partition of unity functions and the denotes 

the number of online eigenvectors that are chosen for each coarse node i. We note that the construction in 







Eq. (14) yields continuous basis functions due to the multiplication of offline eigenvectors with the initial 
(continuous) partition of unity. Next, we define the online space as 


Qon = spanj'i/ii^fc : 1 < i < and 1< k < (15) 

Using a single index notation, we may write Qon = span{'0j}^‘^^, where denotes the total 

number of basis functions in the spaces QoA S'lid is number of coarse mesh nodes. 

Denote the matrix 

i?p = [V’i,...,V'Arf] , 

where ipi are used to denote the nodal values of each basis function defined on the fine grid. 

4.2 Multiscale basis functions for displacements 

For construction of multiscale basis functions for displacements we use similar algorithm that we used for 
pressure. We first construct a snapshot space for each parameter Again, as with pressure we give 

two possible snapshot space choices. 

Snapshot Space 1: As our first possible snapshot space we propose the harmonic extension using . 
We define as the solution to 

L/ cj.snap \ 

a (MUV?/,;- ^,'c) = 0 m w, 

, ( 16 ) 

= ondu:. 

Again, S[^(x) = 6i^k, V/ G and for each fixed parameter values /i^, j = 0,..., N. 

Snapshot Space 2: We could also use the method based on solving an eigenvalue problem with natural 
boundary conditions given by 




CD,snap 


(17) 


Where 


= {C{x]fij)e{ipi),e{(pj))dx, iV„(^j) = / m{x; dx, 

Jn Jn 


and, in the case of linear elasticity m(x; pj) = (A + 2p,e). In a more complicated relation m(x] pj) is related 
to the lower order operators |3]. Again, ipi are the standard fine-scale basis functions, and his is done for 
each fixed parameter values pj, j = 0, 

Define 

U“,p = span{$^7P : 1 < I < k, 0 < j < A}, 

for each coarse subdomain w. We denote the corresponding matrix of snapshot functions, again with similar 
notation, to be 

where Agnap denotes the total number of functions to keep in the snapshot construction. 

Again, we perform a dimension reduction of the space of snapshots by using an auxiliary spectral decom¬ 
position. We solve the parameter-independent eigenvalue problem in the space of snapshots 


where 


= Af A°®$f, 


(18) 


= (^r„ap)^ ^Knap, = (Knap)^ ^K„ap> 

where A and A denote fine scale matrices 


Amn= / {C{x)e{ipm),e{ipn)dx, N^n = / m{x)iprn ' Vn dx. 
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Here, ^pi are fine-scale basis functions. Further, we have 


N N 

c{x) = Y, tjC{x,^j), m{x) = tjm{x] /ij) 
j=i i=i 


is independent of and tj are prescribed non-negative weights. Recall, the main objective is to use the 
offline space to accurately construct a set of multiscale basis functions for each jjLj in the online stage. As 
before for the fluids flow module, at the offline stage of the mechanics the bilinear forms are chosen to be 


parameter-independent, such that there is no need to reconstruct the offline space for each /ij. 

We then choose the smallest eigenvalues from Eq. (181 and form the corresponding eigenvectors in 

the space of snapshots by setting for fc = 1,..., where are the coordinates 


of the vector We denote the span of this reduced space as and denote 


TDU _ 

— 


, _off 


/ 


At the online stage, we use following parameter-dependent eigenvalue problem 




(19) 


where 


A“(m)= f {C{x;n)e{p°P),e{pf))dx={K^fAR:s, 

J to 

= [ rn{x,fi)ip'ppfdx={R^gfNR^g. 

J U) 

Finally, we multiply the linear partition of unity functions by the eigenfunctions in the online space 
to construct the resulting basis functions 

V>^,k = for 1 < z < A^e and 1 < fc < (20) 

where v?™ = denotes the number of online eigenvectors that are chosen for each 

coarse node i. Next, we define the online space as 

Von = span{(/?j,fc : 1 < * < A^c and 1 < A: < N^^}. (21) 

Using a single index notation, we may write Von = span{(/3i}^;^, where = X)£i ^on^ denotes the total 
number of basis functions in the space V^. 

And after construction Von "we denote the matrix 

Ru = [ipi,...,ipN^]'^ , 

where (pi are used to denote the nodal values of each basis function defined on the fine grid. 

4.3 Global coupling 

Now that we have constructed the online spaces for both the fluid and mechanics we now can use this 
parametrized basis at the global level. Indeed, for global coupling we use system of equations to find 
(m”+i\p”+i) e (QonVon), where 

Qon = span{^i}^";^, and Kn = span{(/5j^“;^. 
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Using the matrices 


Rp= [ipi,...,ipNs] : and Ru = , ipN^] , 

we may write matrix analogue for the variational for that will be used for calculation of multiscale solution 
Writing in matrix form, using the notation in (§ , in the online basis we have 

+ aRuGR^pJ^:f+^ = 0 , ( 22 ) 

Rp{S + TB{p^))R^p]Xf+^ + aRpDRluy_^{^+^ = tR^F + R^SR^p^. (23) 


We also note that matrices Rp and may be analogously used in order to project coarse-scale solutions 
onto the fine-grid 


Pltl 


T ms,n-\-l 


T~ti rris'. 

p^pPj+i 


^j+1 


= Ri u 


ms,n+l 

j+i 


5 Numerical Examples 


In this section, we present numerical examples to demonstrate the performance of the GMsFEM for com¬ 
puting the solution of the nonlinear poroelasticity problem in heterogenous domains and complex nonlinear 
dependence on permeability and elastic properties. We use fully coupled scheme for approximation by time 
with Picard iteration to linearize the nonlinearity. We will implement a single complicated geometry with 
contrasting parameter values. Indeed, as noted before, there are many possible nonlinear relations, but here 
we take a an exponential pressure relationship with the permeability. We present the errors with varying 
number of multiscale basis functions and over time for linear and nonlinear case with parameters. 

We proceed as in [9], and we take the computational domain U as a unit square [0,1]^, and set the source 
term f — 0 in §• We utilize heterogeneous coefficients that have different values in two subdomains. We 
denote each region as subdomain 1 and 2, Ui,U 2 , respectively. We use following coefficients: for the Biot 
modulus we take Mi = 1.0, M 2 = 10 in each respective numbered subdomain For permeability we take 
a linear K and nonlinear relation K{p). More specifically, for the linear regime we have 


exp(l) 
exp(10) 


in Ui, 
in U 2 . 


For nonlinear case we use a permeability that depends on pressure p 


K{p) 


exp(p) in Ui, 

exp(lOp) in 172. 


(24) 


(25) 


For fluid-solid coupling constant we have a = 0.9. For the elastic properties we use following coefficients: 
elastic modulus is given by Ei = 10, U 2 = 1 in each respective subdomain 17^, the Poisson’s ratio is ry = 0.22, 
and these can be related to the parameters pi and Ai, for i = 1,2, via the relation 


_ Ei ^ _ Eip 

^^-2{l + r,y (l + ry)(l-2^)’ 

in each subdomain. The subdomains for coefficients shown in Figure where the background media in red 
is the subdomain 1, 17i, and isolated particles and strips in blue are the subdomain 2, 172. 

As we have chosen / = 0 we must use boundary conditions to force flow and mechanics. In these tests, 
as in j9], we use following boundary conditions: 


and 


p = Pi, xGTt, P=Po, x€Tb, = 0, xGFiUFij, 

on 


Ux = 0, ^ =0, a: G F 

dy 


dux 

L, =0, Uy=Q, x€ Fb, 
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Figure 2: Coefficients subdomains. Red is the subdomain 1 and blue is the subdomain 2 


and finally, 


duj; _ 
dx ~ ’ 


9uy 

dy 


a: G Ft U Fij. 


Here F l and Fij are left and right boundaries, Ft and F b are top and bottom boundaries respectively. We 
set po = 0 and pi = 1 to drive the flow, and thus, the mechanics. 



Figure 3: Two coarse grids and fine grid. Left: coarse grid with 36 nodes. Middle: coarse grid with 121 
nodes. Right: fine grid with 3721 nodes. 


In Figure we show the two coarse grids and fine grid. The first coarse grid consists of 36 nodes and 50 
triangle cells, the second coarse grid contains 121 nodes and 200 triangle cells, and the fine mesh consists 
of 3721 nodes and 7200 triangle cells. The number of time steps is Mt = 10 and the maximal time being 
set at Tmax = 0.055. As an initial condition for pressure we use p = po = 0- For the nonlinear solve we use 
Picard iteration for linearization and terminate the iterative loop when \ \pf — PmsIUalf!) ^ ^ = 10“®. 

The reference solution computed by using a standard FEM (linear basis functions for pressure and 
displacements) on the fine grid, Picard type linearization, and using a fully coupled time-splitting scheme. 
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The pressure and the displacement fields on the fine-grid are presented on the left column of Figure and 
Figure 


Fine 


Coarse 



Figure 4: The fine-scale and coarse-scale solutions of the pressure distribution for T = 0.02 and 0.055 (from 
top to bottom) for nonlinear case. The dimension of the fine-scale solution is 11163 and the dimension of 
the coarse space is 864. 


The errors will be measured in relative weighted and relative weighted norm for pressure 

ll„ II _ 

and for displacements, due to the linearity in our Elasticity in this example, we have 


11 ^“ 1 ^ 2 ( 0 ) 


(/f2 ('*'+2 m) (“/.“/Ida:) ^ 
{f[i{'yiuf-Ums),eiuf-Ums))dxy^^ 

{/r!('^(“/).e(«/))da:)^^^ 
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Fine 


Coarse 



Figure 5: The fine-scale and coarse-scale solutions of the displacements and Uy for nonlinear case. The 
dimension of the fine-scale solution is 11163 and the dimension of the coarse space is 864. 


Here {uf,pf) and {ums,Pms) are fine-scale and coarse-scale using GMsFEM solutions, respectively for pres¬ 
sure and displacements. 

In our examples, the nonlinearity resides in the pressure solves. Therefore, we will use the nonlinear 
parameter dependence approach in Section 4.1 For our Elasticity basis construction, we may remain in 


the linear algorithmic approach to construct the online basis. In general, for simulation using GMsFEM we 
first generate a snapshot space using first choice ((10), snapshot space 1) or second choice ((dH), snapshot 
space 2), then we use a spectral decomposition to obtain the offline space, and similarly to obtain the online 
space. For each time step and nonlinear Picard iteration we update the online space for pressure and solve 
the equation (13) utilizing the previously computed solution For construction the snapshot space 2 we 

choose a specified number of eigenfunctions h = 16 for all uji. We select the range of solutions Pmin = 0 and 
Pmax = 1 and divide the domain [pmimPmax] into N equally spaced subdomains to obtain N + 1 discrete 
points po, .■■,Pn- For simulation we use N = 20. 

Recall, we will use a few multiscale basis functions per each coarse node Wi, and these number of coarse 
basis defines the problem size (dimension of online spaces, Qon and I4n)- We suppose that in each patch 
uJi we take the same number of multiscale basis functions for pressure, for all uJi. Similarly 

for displacements we take iV"„ = for all Varying the basis functions in both pressure and 

displacement multiscale spaces we record the errors at the final time. We note that the size of online space 
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and the associated solution accuracy will depend on the number of eigenvectors and that we keep 
in the online space construction. 

We begin first with the purely linear case with K given by (241. In Tables [T] and we present the relative 
weighted and errors for linear case of the coefficients in geometry Figure fusing the fully coupled 
time scheme on two coarse grids. In Table we have a coarse-grid of 36 nodes and in Table we have a 
rehned coarse grid with 121 nodes. We compare these to a fine-scale solution space with dimension 11163. 
In these tables, and are number of multiscale basis functions for each neighborhoods, the second 
column show the dimension of the online space, the next two columns present the relative weighted and 

errors for pressure and last two columns show the relative weighted and errors for displacements. 
We note that as the dimension of the online space increases, because we keep more eigenfunctions Np^, 
in the space construction. We note for the less rehned coarse-grid with 36 nodes the relative weighted 
errors decrease from 36.5% to 0.07% for pressure and from 24.3% to 0.5% for displacements and relative 
weighted errors decrease from 99.0% to 2.7% for pressure and from 37.7% to 3.4% for displacements. We 
note for the rehned coarse-grid with 121 nodes the relative weighted errors decrease from 14.1% to 0.01% 
for pressure and from 26.9% to 0.1% for displacements and relative weighted errors decrease from 82.0% 
to 1.6% for pressure and from 36.1% to 2.5% for displacements. We note that in this example, rehning the 
coarse-grid is not as advantageous to more local basis functions per grid-block. Indeed, with the less rehned 
coarse-grid of 36 nodes and = 12 gives a very good percentage error for a space of dimension 

1296 when compared to the more rehned coarse-grid of 121 nodes and less eigenvectors = 4 with 

space of dimension 1452. 


NP 

on 

dim(Qon; ^n) 

Pressure errors, Sp 
iji 

Displacements errors, 

'on ^ 

2 

360 

0.365 

0.990 

0.243 

0.377 

4 

432 

0.057 

0.435 

0.238 

0.370 

=S 

'on ^ 

2 

648 

0.365 

0.990 

0.108 

0.207 

4 

720 

0.057 

0.435 

0.045 

0.077 

8 

864 

0.001 

0.059 

0.017 

0.072 

= 12 
on 

2 

936 

0.365 

0.990 

0.111 

0.199 

4 

1008 

0.057 

0.435 

0.042 

0.045 

8 

1152 

0.001 

0.059 

0.007 

0.034 

12 

1296 

0.0007 

0.027 

0.005 

0.034 


Table 1: Numerical results for linear problem for coarse mesh with 36 nodes. 


In a similar setting, we consider the nonlinear case of the coefficient with K{p) given by (25). Here we 
will explore the different snapshot spaces available for us in the nonlinear algorithm. Again as in the linear 
case we use two coarse-grids and implement this with a fully coupled time scheme and use Picard iterations 
for the nonlinearity. We present the results in Table for snapshot space 1, the errors are very similar in 
magnitude when compared to the corresponding linear case. In the left side of Table we present the errors 
for 36 nodes in the coarse-grid. The relative weighted errors decrease from 8.1% to 0.09% for pressure 
and from 30.4% to 0.5% for displacements and relative weighted errors decrease from 60.9% to 4.7% for 
pressure and from 38.0% to 3.4% for displacements. In the right side of Tablewe present the errors for 121 
nodes in the coarse-grid. The relative weighted errors decrease from 4.8% to 0.02% for pressure and from 
26.4% to 0.1% for displacements and relative weighted errors decrease from 45.9% to 2.7% for pressure 
and from 35.7% to 2.5% for displacements. For snapshot space 2 we do precisely the same experiment with 
two coarse-grids. We present the errors in Table and again see that the errors are also decrease and have 
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NP 

on 

dini(Qon; ^n) 

Pressure errors, Sp 
iji 

Displacements errors, 
i/i 

N'^ = 4: 

'on ^ 

2 

1210 

0.141 

0.827 

0.269 

0.361 

4 

1452 

0.007 

0.132 

0.240 

0.352 

N'^ =8 

'on ^ 

2 

2178 

0.141 

0.827 

0.069 

0.095 

4 

2420 

0.007 

0.132 

0.024 

0.063 

8 

1904 

0.001 

0.042 

0.015 

0.062 

= 12 

■^'on 

2 

3148 

0.141 

0.827 

0.059 

0.076 

4 

3388 

0.007 

0.132 

0.011 

0.027 

8 

3872 

0.001 

0.042 

0.003 

0.025 

12 

4356 

0.0001 

0.016 

0.001 

0.025 


Table 2: Numerical results for linear problem for coarse mesh with 121 nodes. 


roughly the same behavior. In general, we see that the two snapshot choices in this example do not differ 
greatly and no clear choice arises. In some cases the snapshot space 1 appears to fair better, however, this 
is not always true. Finally, we note that, for solution of nonlinear problem in each time step, the Picard 
iteration converges after about 3 steps. 

To show the stability of the multiscale spaces over time we include time plots. We include plots over 
time of the error with respect to number of basis functions used. To get an idea of the behavior we only 
present the results for snapshot space 1 for two coarse grids. In Figure and we show errors over time for 
Non = -^on = N^n = 4, 8,12, respectively. We observe that errors decrease as we increase the dimension of 
the offline space as expected and the basis appears to be robust with respect to longer times. 


NP 

on 

a 

P 



N'^ = 4 

'on ^ 

2 

0.081 

0.609 

0.304 

0.380 

4 

0.019 

0.242 

0.254 

0.371 

=8 
'on ^ 

2 

0.082 

0.607 

0.091 

0.104 

4 

0.021 

0.241 

0.023 

0.074 

8 

0.001 

0.087 

0.016 

0.072 

= 12 

■^'on 

2 

0.082 

0.607 

0.085 

0.077 

4 

0.021 

0.241 

0.013 

0.037 

8 

0.001 

0.087 

0.007 

0.034 

12 

0.0009 

0.047 

0.005 

0.034 


NP 

on 

£ 

P 



=4 

'on ^ 

2 

0.048 

0.459 

0.264 

0.357 

4 

0.008 

0.132 

0.235 

0.351 

N'^ =8 

'on 

2 

0.048 

0.457 

0.063 

0.079 

4 

0.006 

0.130 

0.022 

0.063 

8 

0.001 

0.053 

0.015 

0.062 

N'^ = i: 

-^'on 

> 

2 

0.048 

0.457 

0.052 

0.051 

4 

0.006 

0.130 

0.009 

0.026 

8 

0.001 

0.053 

0.002 

0.025 

12 

0.0002 

0.027 

0.001 

0.025 


Table 3: Numerical results for nonlinear problem using snapshot space 1. Left: for coarse mesh with 36 
nodes. Right: for coarse mesh with 121 nodes. 
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Figure 6: Weighted are on the top and are on the bottom. Errors for pressure are on the left and 
displacements are on the right for nonlinear problem on coarse mesh with 36 nodes. 


6 Conclusion 

Modeling and simulation of a nonlinear poroelastic media is challenging due the heterogeneities and the 
nonlinear dependence on the coefficients. Thus, in this paper we developed a Generalized Multiscale Finite 
Element Method for a nonlinear poroelastic media. We gave a general nonlinear poroelasticity model in 
the framework of the Biot equations, where we had possibly complex nonlinear dependence on permeability 
fields and elasticity tensors. As the Nonlinear GMsEEMs treat nonlinearities as a parameter, we linearize 
the equations in a time-staggered Picard iteration formulation. We then outlined the construction of the 
multiscale spaces offline and online spaces. The algorithm is then implemented on a single geometry with 
two different cases of permeability fields. The first being the standard linear case and a second nonlinear 
relation depending on pressure where a parameter spaces are considered with offline and online spaces. We 
presented the errors relative to the fine scale solution with varying multiscale basis functions and coarse-grid 
refinements. Einally, we showed the robustness of the modes for longer time simulations. 
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on 


Sp 





^on = 

1 


2 

0.063 

0.551 

0.302 

0.380 

4 
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0.256 

0.267 

0.372 



N'^ = J 

'^on ^ 
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